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AN ANALYTIC METHOD TO ACCOUNT FOR DP \G 
IN THE VINTI SATELLITE THEORY 

J. S. Watson 
G. D. Mistretta 
N. L. Bonavito 

ABSTRACT 

In order to retain separability in the Vinti theory of earth satellite motion 
when a non-conservative force such as air drag is considered, a set of varia- 
tional equations for the orbital elements are introduced, and expressed as func- 
tions of the transverse, radial, and normal components of the nonconservative 
forces acting on the system. In this approach, the Hamiltonian is preserved in 
form, and remains the total energy, but the initial or boundary conditions and 
hence the Jacobi constants of the motion advance with time through the varia- 
tional equations. In particular, the atmospheric density profile is written as a 
'fitted 1 exponential function of the eccentric anomaly, which adheres to tabular 
data at all altitudes and simultaneously reduces the variational equations to 
indefinite integrals with closed form evaluations, whose limits are in terms of 
the eccentric anomaly. The values of the limits for any arbitrary time interval 
are obtained from the Vinti program. 

Results of this technique for the case of the intense air drag satellites San 
Marco-2 and Air Force Cannonball are given. These results indicate that the 
satellite ephemerides produced by this theory in conjunction with the Vinti pro- 
gram are of very high accuracy. In addition, since the program is entirely 
analytic, several months of ephemerides can be obtained within a few seconds of 
computer time. 
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INTRODUCTION 

Vinti (Reference 1) has shown that if a satellite orbit is described by meanE ^ 

of osculating Jacobi a’s and 5’s of a separable problem, then a perturbing ■ 

i 

force F makes them vary according to J 

i 

1 

? 

= F • 37/3/^, - - p • dT/dc^, (K - 1, 2. 3)- j 

! 

Here r is the position vector of the satellite and F is any perturbing force, ! 

„ j 

conservative or non-conservative. If F is the force of air drag, the interaction ; 

of drag with oblateness makes it desirable to obtain variations to the order 
drag x J 2 , where J 2 is the coefficient of the second zonal harmonic of the 
Earth’s gravitational potential. The physical reason for carrying these deriva- 
tives through order J 2 is the strong variation of drag with perigee height. In 
the present paper, we have been able to account for this effect without intro- 
ducing the J 2 terms into <57/?/3 K and 37/3a R . The logic behind our approach 
requires a rather careful exposition which we shall go into in detail in Section I. 

The essence of the method is that for a given time interval, one always does 

both a drag free calculation, and an oblateness free calculation with drag, and 

that these two calculations are done in a self-consistent iterative manner such 

that the mean orbital elements never go far astray. The appropriate criterion, 

to make sure that the drag-oblateness interaction is being properly accounted 

for, is that the perigee height corresponding to initial and final orbital elements : 

of a given interval shall not change by more than some predetermined amount. 

The complexity of those papers which attempt to handle the oblateness- 
drag interaction in a straight forward manner (References 2, 3) illustrates the 
desirability of finding a new approach. That is the purpose of this paper. i 


L 
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I. STATEMENT OF THE PROBLEM 

In this paper we consider the motion oi an artificial Earth satellite in the 
presence of air drag and the Earth's gravitational potential. In contrast to the 
classical methods of numerical integration, our approach will be to present a 
quadrature algorithm employing analytical expressions for the variation of 
orbital elements produced by air drag. These expressions are well-defined 
over expanded subintervals of the solution, and produce accurate agreement 
with profiles of tabular density. This procedure then allows a flexibility in the 
selection of end points of the subintervals, which in turn insures a minimum 
error bound on the required analytical function. For convenience we shall 
henceforth refer to the algorithm as the BMW (Bonavito-Mistretta-Watson) 
aerodynamic method. In this method the effect of oblateness is accounted for 
by the Vinli Spheroidal Theory (Reference 4). The changes due to atmospheric 
resistance for a non-rotating sphere are accounted for by the solutions of the 
variational equations without oblateness (Reference 5). 

Normally, one would wish to represent the variation of atmospheric density 
by an exponential whose power is a function of the difference between the satel- 
lite height and the altitude at a predetermined density (Reference 6). Such a 
representation is usually valid only in a neighborhood of this boundary value. 

The neighborhood or region over which this density representation is in agree- 
ment with tabular data such as provided by the U. S. Standard Atmosphere Sup- 
plements, 1966, (Reference 7), is one in which the density scale height is observed 
to vary in an approximate linear fashion. Throughout our calculations a set of 
such regions is chosen to meet this requirement. In addition, the initial or 
boundary value of the atmospheric density for each of these regions is also 
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supplied by the Supplements. Such a representation for atmospheric density 
would, except under special conditions, exclude closed form integration of the 
variational equations. To avoid this difficulty, we approximate the atmospheric 
density variation by an expression which is made to adhere closely to the 
numerical values of the aforementioned model. By adjusting or advancing 
boundary conditions over several selected arcs or layers of atmospheric density, 
we produce a profile that closely agrees with the tabulated data for all heights. 
The degree to which our results compare with tabular values (Spring- Fall model, 
1100°K exospheric temperature) from the U. S. Standard Atmospheric Supple- 
ments, 1.966, is shown in Table I. Between the heights of 205 kilometers and 650 
kilometers, the discrepancy is less than two percent. From symmetry consider- 
ations, these sets of boundary conditions for one representation can be deter- 
mined during the first half revolution of an orbit. These density variation pro- 
files then are held fixed until such time as the perigee height changes by some 
predetermined amount. At this point the hours dar.y-eonditions are redetermined 
over the first half revolution away from perigee, corresponding to a chosen 
epoch. This again produces a total density profile that is in close agreement 
with the values from the Supplement tables. Experience indicates that given a 
criterion of one kilometer change in perigee, this redetermination is not neces- 
sary for nearly two months in the cases of the San Marco-2 and Cannonball 
satellites, but becomes more frequent near the end of the lifetime of each 
spacecraft. 

In our final expressions for these variations, rotation of the atmosphere 
Is accounted for while an oblate atmosphere is considered in a partial fashion 
by including the oblateness narameter in the exponential term of the density 


i 

K 

f 

l 

I 

i 

( 

} 




variation. Although the variational equations do not contain the effects of oblate- 
ness explicitly, the interaction between oblateness and drag is accounted for 
implicitly during the computational procedure in the following manner (See 
Figure 5). 

1. From the initial conditions calculate the Izsak elements of the Vinti 
Spheroidal Theory corresponding to that eopcli. 

2. With these, together with tabular data or. air density, obtain the density 
variation profiles and corresponding values for changes in the elements arising 
from air drag for those altitudes during the first half revolution past perigee 
and beyond epoch, corresponding to these density profiles. 

3. For any desired time interval, calculate using oblateness only (without 
drag) 7 and 7 and the value of the eccentric anomaly E u at the end of this 
interval from the initial given set of Vinti or Izsak elements. 

4. For this same time interval, calculate from drag only (without oblateness) 
the total change in the original or given set of Vinti elements, using the value of 
the eccentric anomaly obtained in step (3) as the upper limit in the analytical 
expressions for these changes in the elements arising from drag. The total 
change in the elements from epoch to this time, produced by drag only, are 
obtained by multiplying the revolutions to this point in the orbit by the sum of 
the individual (analytically expressed) corrections or changes to the elements 
obtained in the first half revolution. 

5. Add these changes arising from drag of step (4) to the epoch values of 

the Vinti elements. . . 

6. For the same time interval, repeat step (3) with the new elements ob- 


tained in step (5). 
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7. Compare the last two values of eccentric anomaly. If their difference is less 
than some arbitrary preassigned small positive number e , then the differences 
between the original noniterated set of elements ai°> , /?(°) (i = 1,2,3), and the 
iterated set a< n >/3f n > approaches a constant. That is, when 


then, 


£("> _ E( n- n 




< £. 


&a = (a f0 i - — constant . 

Similarly, 

A/3 = (/3< 0) - /5 (n> ) = G(/2 (n-s> , E< n ' ] >') - constant 

where (n) is the iteration number. Note that the functional representation of 
the variational equations solution on the right hand side contains the previous 
iterated values of the elements and the eccentric anomaly, (a (n-1 l , /3< n-1 l, 
E <n-n ^ 

8. If the above criterion on the eccentric anomaly upper limit is met, then 
we accept the a< n l , as the new Vinti elements that describe the orbit from 
the original epoch to the end of the given time interval. 

9. If the criterion is not met, repeat steps (3) through (7) always utilizing 
the iterated values of the Vinti elements in the ea^ulation of step (3) until such 


F( n > _ ir( n "i) 

U U 

£("> 

u 


time as 
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In functional terms, if by x we denote the elements of the Vinti Satellite 
Theory and f(x, t) represents the Vinti Spheroidal Method (oblateness only) 
solution which for a given time interval yields a value of the eccentric anomaly 
E u corresponding to a time t at the end of that interval, and g[x, f(X,t)] repre- 
sents the correction arising from drag plus the oblateness calculation; then for 
a fixed value of time t, the above algorithm is an iterative solution for x using 
the equation 

f(x, t) = f(g.fx, f (x. t)l, t) 


in a self-consistent method. Thus, for a fixed time interval 
- E<” - (Aa”>. AA”>> - ... - (a'"' 11 - 


- E (n) - (Aa<"\ A/3 (n) ) - (a ( "\ /3 <n) ). 

U drag 


This algorithm converges when 


t.<") 


< € . 
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II. AERODYNAMIC DRAG VARIATIONAL EQUATIONS 

From Sterne (Reference 5) , the equations for variations of the orbital ele- 
ments due to air drag and without oblateness are given by 


da 

dt 


/ a 3\l/2 1 

2 J o ' 2 [ Re sin v + T (1 + e cos v)j 


de 

dt 

di 

dt 


dF 


dfl 

dt 


dM 

dt 


/ a * 

(7I ) [R s i n v + T( cos v + cos E)] 

( r cos <b \ 


(i) 


tjz 5 « (i\ - e 2 (l + — J sinv\ , 

r. 1 *1 ~ e cos v l , .J VP/ \ f r cos i sin \p \ 

^ "*• 1 \ — « y ■*(„.> vi-.4 ,i„i) 


/ r sin </* \ 

\na 2 t'i - e 2 sin i / 

/(l - e 2 ) 

n + R l nae cos v ' 


2r 


- T 


r {l - e 2 )(l + sin v V 


nae 


Here \p is the argument of latitude, r = a(l-ecosE), p = atl-e 2 ), and 
n - fi 1/2 a _3/2 , where pt. is the product of the gravitational constant and the sum 
of the masses in the two body problem. E and v are the eccentric and true 
anomaly respectively and I is the inclination. The drag perturbing force is 
resolved into the following components: 

R is in the direction of the position vector from the force center to the satellite. 
T is perpendicular to R, lies in the orbital plane, and is positive in the direc- 
tion of motion. 



10 


W is mutually perpendicular to R and T and completes a right handed set of 
component directions . 

In terms of the eccentric anomaly, these are given by 


R = - - C D — aepV sin E 

2 B m dt 

T = - Ir A (i - e 2 ) 1/2 a PV 
2 u m 



(1 - e cos E) 2 

(1 - e 2 ) - 


dE 

dt 


( 2 ) 


W = - - CL - apa> V/J.~ 1/2 £.3/2 
2 m * 


x (1 - e cos E> 2 V*S cos \p 


dE 

dt 


Here, is the angular velocity of rotation of the Earth, P is the atmospheric 
density. A is the projected area of the satellite, ra is the mass of the satellite, 
and C D is a drag parameter. S is a parameter related to the orbital inclina- 
tion and is given by 

S = sin 2 i. (3) 

The velocity of the satellite relative to the atmosphere can be written in vector 
notation as 

V = T - S x V, 

Neglecting the term of order co * in the magnitude of V to be employed in (2), 
one obtains the expression given by Sterne to be 

V = V /2 (1 + e cos E)" 1/2 (1 - e cos E)‘ 1/2 [(1 + e cos E) - d(l + e cos E)] 


{ 

I 

i 
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where 


d = ai g) u. 1/2 a 3/2 (l - e 2 ) ,/2 (1 - S) l// sgn a 


1/2 


where sgn a = +1 or -1 according to whether the orbit is direct or retrograde. 
From (3) above we have 


ds „ . di 

=2sinicosi 

dt dt 


or 


f = 2 YsYT-~sf t 




§ 


§. 


& 


Inserting di/dt from (1), we can then write the variation of the inclination 
parameter as 


dS _ 2 yS(l - S) rw COS l p 

Let us now write the right ascension of the ascending node and the argument of 
perigee as J3^ and respectively. In addition, the mean anomaly is related 
to the time of perigee passage/^ by M = n(t + p±). Differentiating, we have that 

d/3 i _ 1 fdM 


where 


d/3 i _ 1 fdM ft x B \ dn nl 

-dF'n [d^~ ( J 


^ _ 2.u 1/2 a' s/2 If 

dt 2 dt ’ 


Using the above results, together with (2), the solutions of the variational equa- 
tions due to air drag without oblateness are given from (1) by 


-3- -i 


a ■ 


I 
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= C 

m t 


p{(l - e cos E) 1/1 (1 + e cos E)~ ,/2 [(1 + e cos E) 


- d(l - e cos E)] 2 } dE 


Ae = - B (1 - e 2 ) f E p (i± ° cos E V V2 [l - d < 1 - 

m J \1 - c cos E/ L (1 + 


e cos E) 
e cos E) 


x jcos E - ^ (1 - e 2 ) 1 (1 - e cos E) (1 cos E - e - e 2 cos 2 e| dE 


AS = - - - D - w — (1 - e 2 ) 1/2 S Vl - S f p(l -i e cos E) 1 2 
m n I 


x (1 - e cos E) S/2 cos 2 </<[(l + e cos) - d(l - e cos E)j dE 


A^ 3 = C D Ap" l/2 a s/2 o) s (l - e 2 )~ 1/2 J 2 p sin cos </> 


x (1 + e cos E) 1/2 (1 - e cos E) 2 [(1 + c cos E) - d(l - e cos E)] dE 


a ^ 2 = - ~c D -ae _, (i - eV i/2 

2 m 


r h- 


e sin 'p cos i p — 1 yi - S sgn a. 


x (1 + e cos E) 1/2 (1 - e cos E) 5/2 [(1 + e cos E) - d(l - e cos E)] 


+ 2(1 + e cos E) 1/2 (1 - e cos E) 1/2 sin E[(l + e cos E) - d(l - e cos E>] 


x Jl - e2 - — (1 - c cos E) (2 - e 2 - e cos E)J > dE 
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and 


40, = J tl (t + [i ) - _ <A0 a + /TTS A/3 3 sgn a 3 ) 

2 a 1 n 


C.A 


?r 


p( 1 - e cos E) 1/2 (1 + e cos E) 


- 1/2 


x f(l + e cos E) - d(l - e cos E)] sin EdE. 

Here Ej and E 2 are the values of the eccentric anomaly at arbitrary times tj 
and t 2 respectively and the time t is usually taken to be the time of epoch. 

The solution of (4) will be considered in Section IV after describing the form 
of p in Section ill. 


HI. ATMOSPHERIC DENSITY REPRESENTATION 

At this point we consider the expression for atmospheric density variation 
given by King-Hele (Reference 6). An expression for the atmospheric density 
as a function of the eccentric anomaly is given by 

p = p p {l + b(r - r p ) 2 } exp f- — - - ? 


(1 - cos E)| (5) 

where x = ae, H p = density scale height at perigee which in equation (5) is 
assumed to vary linearly with the perigee height. 

Let us now write equation (5) in a more general form so as to represent 
the atmospheric density variation over any interval within the orbit in which 


= p p { 1 + bx 2 (l - cos E) } exp 


R 
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the boundary conditions are known. A generalization of (5) is desirable for our 
purposes to insure that the resultant analytic atmospheric model rigidly adheres 
to a tabular set of densities at all altitudes. To obtain this generalization, let 
the subscripts L and U refer to a lower and upper point of an interval respec- 
tively. This interval [ E L , E y ] is judiciously chosen (by a method discussed 
later in this section) and allows either 

0 < < 7T and 0 < < rr 


or 


7T < < 2n and n < < 2rr. 


Letting 


and 


r L = a (1 - e cos E^), = a(l - e cos Ey). 


r rain = min < r L' r u>* r ma* = maX < r L’ r„), 
permits p to be given by an expression of the type of equation (5) over the re- 
stricted domain [El, E,j] by 


p(r) = p L { 1 + b(L, U) (r - r mjn ) 2 } expf- 


^ ! r < r < r 


^ J 


or 


p(E) = p L {l 4 b(L, U) x 2 (cos Ej^- cos E) } exp|- — (cos - cos E)^. El<E<E u ( 6 ) 


’K 


>} 


where p L and Hl are available from density tables, El and E^ are available 
from the orbit theory (as are r L and r^, and b(L,U) is derivable from (6) by 
setting pCEy) = p u which yields b(L, U) = b(L, U; p L . P v - r L” 
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If equation (6) above is substituted into (4) and the resulting expressions 
are expanded in powers of the eccentricity e, this yields terms of the form 
r E 2 

I s i n n E cos" E exp{- (x/H^) (cos Ej^ - cos E)} dE 

Je, 

where m 1 0, nS 0, and H L is held constant over the range of a subinterval of 
integration and x is held constant over the entire range of integration. For the 
case n= 0, this expression can be evaluated in terms of Bessel functions for the 
limits (0 ,tt). The parameter n, however, does assume the value zero; and, in 
addition, it is desirable to alter the form of the integrand so as to insure that 
the expression under the integral sign is integrable over any interval while 
simultaneously maintaining its theoretical physical content. In short, the vari- 
ational equations cannot be solved by th» Fundamental Theorem of Calculus 
since an antiderivative of sin" E cos m E exp { ~(x/ll L ) (cos - cos E)> , (n ^ 1), 
is not expressible in terms of elementary functions. 

To achieve this purpose let us now consider the following: rewrite the 
expression for atmospheric density (5) as 

p(E) = p L [1 + b(L, U) x 2 (cos Ej^ - cos E) 2 ] h t (E) exp[b 2 (E) -E] 

El - E - Ep (7 > 

where b, (E), b 2 (E) are arbitrary functions of E which will be chosen to pre- 
serve closed form integrability of the variational equations (4) while retaining 
near complete numerical agreement with the exponential term of King-Hele's 
expression for p (E). 

The selection of the two functions bj(E), b 2 (E) between an upper limit of 
integration E 2 and a lower limit E } will be made in the following manner. 


/ 



I 
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Consider the interval [E*, E£] subdivided into a collection of n nonoverlapping 
subintervals I, , (•£ = 1, 2, . . . , n), defined by the partitioning 

e; = x 0 < x, < < x n _, < x n = e;. 

For any subinterval = [x,£ , x^] , define 


for 


*£-! i E £ X; 


b,(E) = 
b, (E) = d, 

where c*, d^ are constants. The determination of cy , will be accomplished 
by applying the technique of classical weighted least squares. In this manner 
the definition of b,(E) and b 2 (E) over the domain [E*, E* ] are the step functions 


ME) = c, x Q < E < Xj 

= c 2 Xj < E < x 2 

= C n X r-1 ^ E ^ X „ 

b 2 (E) = dj x 0 <E<x, 

= d 2 x l < E < x 2 


ad x , < E < 

n n- 1 •“ — 
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Consider observational data Y t generated within by repeated evaluation 
of the function 

h(E. ) = exp[-&(cos Ej_ - cos E £ )] 

where Cl = x/Hj^ and x£_j <E_< xj>. To define the pseudo-regression situation 
governing our estimation problem, consider the implicitly linear, exponential 
model 

Y. =fl 1 e. expCd 2 E,) i = 1, 2, . . . , ro < 8 > 

relating the independent variable E and the dependent or response variable Y 
that, unlike the general regression response variable, displays no random vari- 
ability. Hence, the random variables e. denoting the dispersion characteristics 
of Y. about the theoretical regression line become meaningless in the curve 
fitting analysis. That is, the random variables e. are degenerate in the sense 
that their probability density functions pfCj ) are given by 

P( £ i)=l £ i =B i 

= 0 otherwise 

where B. can be considered a "fitting bias" at E = E. since 

E(€ i ) = B i Var(c. ) = 0. 

While all practical statistical properties of the regression analysis become 
lost, the technique is not degraded as a numerical tool for approximating with 
great precision complex functional forms over restricted regions with relatively 
simple fractions. 





■ i 1111111 imi .uni ii 
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To display the implicit linear form of (8), we take the natural logarithm of 

both sides of (8) and obtain the equivalent form 

lnY. = ln/8 1 +/3 2 E i + lne. 
or 

Y! =/S; +/5 2 E. + e! (i = 1, 2, ---..n). 

If nonnegative weights a> 1 , f ... , <^ m are available and y? are evaluated 

from y! = In h (E. ) = - Q (cos x^_j - cos E. ), then the well known result 




I 

are the values of ft lt fi 2 that minimize the v, 'sighted sum of squares 

m 

i i 

Thus, bj = exp (b’) and b 2 are the classical weighted least squares estimated 
of , /3 2 respectively. 

The selection of the mesh points over an interval of integration [E* , E*] 
is obtained through an automated search approach utilizing the weighted least 
squares process previously described. This technique selects subintervals of 
maximum size while retaining a user selected error tolerance between the true 
and predicted function. ■ 



Furthermore, if one sets = 0, EJ = v, and derives a 6et of coefficients 
for the n selected subintervals in the manner previously defined, the functions 



i 
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h(E) and p (E) will have been fitted for all values of E over which x has been 
assumed constant. This is true since h(E) is a periodic function with period 2n 
and in the fundamental period (0, 2 ti) is symmetrically distributed about the 
axis of symmetry E = v . 

With the use of more general expression for density (p), the variational 
equations (4) integrated between E’J and EJ where 0 < E* < EJ < 77 » take on the 
following general functional form 

n- 1 - x 4 

A q k = £ J , + 1 it j + I exp(d j4l E)} dG k (E) (k = 1, 2, •••.6) 

jnO 

where dG k (E) = g k (sinE, cosE) dE. Having fitted h(E) from (0 -rr), the inte- 
gration of the variational equations between any two arbitrary limits E*> 0, 

Ejf > E 1 }, lakes the general functional form 


n* 1 /*j( t 

Aq k = m j + i J J+ {t j + i ex P Cd j + 1 E > dG k <W 

•ist» * i 


n-^ /* 2 w_x 

2_. m 2n-i J ' {t 2„-j eX P [d in-j E] } dG k< E > 

jTo 277 " x i + i 


CE*niod in 

J t* exp (d£ 


1 E) dG k (E) 


( 9 ) 




t* exp (d* E) dG k (E) (k = 1, 2. •••, 6) 


Ej mod in 


t 




where , ft = 1, 2, . . . , 2n), are nonnegative integers; tj , t k ', d k (j = 1, 2, . . . 
n), (k = n + 1, n + 2, . . . , 2n) are functions of c ■ , d . ; f is an integer between 0 
and 2n - 2; t; is an integer between 1 and 2n - 1; the grid points X n+J , x nt2 , . . . 
x 2r are defined from x 0 , x t , . . . , x n ; and (t*, t*) may be either (t^, t*), 

(t£, t„), (t F > tl,), (t^ , tj) depending upon the values Ef mod 277 , E* mod 2 -,t, The 
logical structure of (13) and the determination of the above parameters is too 
lengthy to be given here but is presented in Appendix for completeness. It must 
be emphasized that (13) is general for all computations, but is valid only while 
x is assumed constant, whereby a new fit is obtained and new constants are 
determined to integrate (4) by the general form (9). 

Using these results we can now solve the density expression for the value 
of b within each of the selected or fitted intervals from the expression 

-P, + exp[(r - r )/Hj 
b(L, U) = — — - — 

^L< r u " r L> 2 

where the subscripts U and L still refer to the upper and lower points of an 
interval. For example, p is that value of density at the initial point of the 
interval which is known from the tables, and the corresponding density scale 
height H l can then be computed. The end point r 0 is known as a function of E, 
Therefore 1^ is also ’advanced* in a similar fashion as p L during the fitting 
of the exponential. The above method of fitting King-Hele’s expression over 
several intervals of an orbit is sufficient to ’reproduce’ the tabular model 
atmosphere with a discrepancy of less than two percent, and simultaneously 
yield closed form integrable equations for the variation of the elements. 
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Table I provides a compilation of an estimated versus a standard model 
atmosphere. The estimated densities eiven in column 3 were obtained via equa- 
tion (7). The standard model atmosphere given in column 2 was obtained by 
using an algorithm which generated as closely as possible specific tables given 
in Reference 7. The standard model atmosphere used here closely resembles 
the Spring- Tall model with an Exospheric Temperature of 1100° Kelvin. The 
first column of the table is given in kilometers, and the remaining three col- 
umns are given in grams per cubic kilometer. 

IV. SOLUTION OF VARIATIONAL EQUATIONS 

In light of the analysis done in Section m above, we now return to the vari- 
ational equations (4). If we now combine the radical terms in (4), we obtain a 
set of expressions containing the forms (1 ±x) n , and (1 ± x) ", where n = 1/2 
and x = e 2 cos 2 E. If one assumes that x will never get too large, that is, for 
drag satellites, e will not be larger than 0.2, then the above terms can be ex- 
panded as 

(1 ± X T = 1 ± ny 4 n - (n 2 ~ X 2 . 

and 

(1 ± x)‘ n = 1 7 n X + n(n 2 * X 2 - < 10 ) 

Adopting expression (6) for atmospheric density variation, the variational equa- 
tions are then reduced to solving a set of indefinite integrals of the form 

f E 2 s 

kj I exp{-(x/R L ) (cos 1^ - cos E)} cos E sin p EdE 

Jr.. 


( 11 ) 
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where kj is simply the constant coefficient of each corresponding integral. If, 
on the other hand, the eccentricity is somewliat Larger than normal, say around 
0.5 or so, then it will be necessary to include several more terms in the expan- 
sion (10) above. This however is an algebraic problem, and computationally 
speaking has only the effect of changing the overall coefficients, { k. } , in the 
final algorithm. The integrals are still of the form (11) above. Utilizing the 
fitting described in section U, the exponential of the density can now be repre- 
sented as 

b,(E) exp [b 2 (E) ■ E] 

in which b 2 (E) and b 2 (E) are determined for one or several segments within an 
orbital revolution. Results indicate that b } (E) and b 2 (E) remain fixed for a 
considerable length of the satellite's lifetime. Refitting becomes necessary only 
when a and e change appreciably. This is found to occur more frequently near 
the end of the satellite's lifetime. In any event, the fitting procedure is instan- 
taneous, and the calculation proceeds without interruption. 

If in the drag acceleration one desires to express the density exponential 
so as to involve the product of J 2 , the coefficient of the Earth's second zonal 
harmonic and the air density, one can follow the procedure of Sherrill (Refer- 
ence 3) and write this term as 

exp J-(x/H ) [fl - e cos E) + — ( 1 — ■ T) — -) + ccr -q 7 

\ L 2a 2 V 1 - e cos E/ E . 

where c 2 = J 2 R 2 . Here R e is the value of the Earth's equatorial radius, <r E is 
the equatorial radius of the oblate spheroid which passes through the initial 
perigee point of the satellite, and the flattening e = 1/298. Since this addition 


1 



23 


f 

■ < 
I 
r 


E 

r 


fir 




f 

££ 

r. 


£•; 

jfej 

& 


»■ 

& 

#?-■ 


is accounted for rather simply by the fitting procedure, the cross term is re- 
tained in the exponential. Using these results, equation (li) is rewritten as 

" E 2 


V E > J 


exp[b 2 (E)'E] cos' E - sin p EdE. 


( 12 ) 


Here p takes the values 0, 1, and 2. When p = 0, 4 takes on all values from 
zero tluough twelve. When p = 1, goes from zero through eleven; and when 
p = 2, ranges from zero through ten. For the case p = 0, expression (12) re- 
duces to 


f r 

[q = b,(E) I exp[b 2 (E)-E] cos 1 ' E dE . 


When -t is even ( -J, = 2n) , this becomes 


r2n . 


bj(E) 


2?n 


exp[b 2 (E)*E] 


(n ! ) 2 b 2 (E) 


(13) 


n- 1 


2.xp M E,.E)^ (2WTJTTP- 


i b 2 (E) cos (2n - 2 j ) E + (2n - 2 j ) s in(2h - 2 j ) E 


i ■«» 


b 2 2 (E) + (2n - 2j) 2 


When •£ is odd (-£ = 2n+l) we have 


b. (E) 

T 2n+1 _ lV / 

A o 

2" 


j -0 


[b 2 (E) cos (2n - 2j + 1) E + (2n - 2j + 1) sin(2n - 2j + 1) E 
bf(E) + (2n - 2j 7 l) 2 


( 14 ) 
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For p = 2 expression (12) becomes 


^2 “ *i( £ > 


p2 

JL, exp 


[b,(E)-E] cos^ E sin 2 EdE 


j exp[b 2 (E) - E] cos^EdE- 

r z 2 

1 exp[b 2 (E)-l 

I] cos^ + 2 EdE 

Jv 
L i 




The calculation here is identical to that for I % with results given by equations 
(13) and (14) for 4 even and odd. 

Finally, for p = 1, we have from (12) 

s 


r E 2 


1^ = bj(E) I exp[b 2 (E)-E] cos' Esin EdE. 


An integration by parts reduces this to 


i; = bj(E) 


-exp [b,(E) • E] cos^ + I F. 


•t + 1 


Sf 


2 o 

exp[b,(E)-E] cos +1 EdE 


Here again, the integral part of Ij is given by equations (13) and (14) for the cases 
that (i + 1) is even or odd. Let us now define the following: 

k 0 = bj (E) p L (l + b(L, U) a 2 e 2 cos 2 E^ 

*1 ~ ~ 2k 2 cos E l 

= b t (E) p L b(L, U) a 2 e 2 

c{ = 0. C’ = 0, C* = 1, C* = 2e, 


C 1 - - e ? C 1 - e 3 r 1 = — e 4 C 1 - - e s 

'“S o ’ U 6 e ’ ^7 o e> '“8 - T e ’ 
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Using the above quantities, we now write the solutions of the variational equa- 
tion (4) as 
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1 C n Aa 

Ae = - i r -2(l - c 2 ) I, - d(2(l - c 2 ) I ]0 

+ 2c” 1 j( g - e(I 4 T I ? )) 

+ d 2 (2e“ 1 I 1(j - e(I 5 + I„))]. 

CL A , 

AS = a 5/7 fi~ l/2 (1 - e 2 )" 1/2 S yi - s sgn a, 

m 

{cos 2 /3 2 (I 7 - 2I 8 ) + e 2 cos 2 fi 2 I 4 

- 2(1 - e 2 ) 1/2 s in [i 2 cos /Sj (I 2 - el, 2 ) 

+ (1 - e 2 ) sin 2 f>> 2 1 13 - d(cos 2 /i 2 (I, 4 - 21, s ) 

+ e 2 cos 2 fij I s - 2(1 - e 2 ) ,/2 sin cos /? 2 ( 1 1 7 
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V. RESULTS 


In order to test the BMW theory, two heavy air drag satellites with some- 
what dissimilar characteristics are considered. These are the Italian San 
Marco-2 and the U.S. Air Force Cannonball (OAR- 901). Data on these spacecrafts 
is as follows: 

San Marco- 2: 

Mass m = 129.27383 kgm. 

Projected Area A = 3425.3397 cm 2 . 

Drag Coefficient C D = 2.1 
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Initial epoch: April 26, 1967, 10 hours, 12 minutes. 

Initial conditions (Inertial Cartesian): 

x = +0.58725272, x = -0.82890608 

y = +0.84923499, y = -rO.56396878 

z = -0.05068537, z = +0.01219124 

Here, x, y, and z are in units of earth radii (6378.166 kilometers), and x, y, and 

z in units of earth radii per canonical unit of time (806.812 seconds). 

Apogee height = 736.00 kilometers 
Perigee height = 205.60 kilometers 
Eccentricity = 0.0387 
Inclination = 2.87 degrees. 

Cannonball: 

Mass m =■ 362.87392 kgm. 

Projected Area A = 3423.6195 cm 2 . 

Drag Coefficient C u = 2.1 

Initial Epoch: August 7, 1971, 0 hours, 20 minutes. 

Initial Conditions (Inertial Cartesian): 

X = -0.9428339, x = -0.26977565 

y = -0.28161629, y = -0.40945775 

z = +0.28038380, z = -1.0092100 

Apogee height = 5 957.20 kilometers 
Perigee height = 130.16 kilometers 
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Eccentricity = 0.1230 
Inclination = 92.00 degrees. 

In both cases, the boundary values of the atmospheric density profile 
given by p L are taken from a static model, namely, the 1966 U. S. Standard 
Air Force Supplements. The profiles used here include a Spring- Fall 
model with an 1100 degree exospheric temperature, a Winter, 800 degree 
model, and a Summer 1000 degree model. While these profiles are 
adequate if chosen carefully, it is felt that a somewhat mors sophisticated and 
dynamic model such as given by Jacchia (Reference 8) would not only improve 
results but also make them more reliable. 

Figure 1, shows the variation over one orbital period (from time of in- 
sertion) of the semi-major axis, eccentricity, and (z 3 2 + £ 3 ) respectively for 
San Marco -2. Here, a Q and e Q are taken to be the initial values of these 
parameters. For the semi major axis, we have initially a secular decrease of 
108 parts in 10 6 per revolution, with a periodic variation superimposed, of about 
25 parts in 10 6 with the orbital period. The eccentricity shows a secular de- 
crease of approximately 97 parts in 10 6 per revolution, and a periodic variation 
superimposed on it with an amplitude of about 22 x 10" 6 , with the orbital period. 
Figure 1 shows that (yS 2 + fi 3 ) shifts back and forth hy - out 8.3 seconds. Figure 2 
is a similar graph for Cannonball. It is seen that the variations for a and e 
here are somewhat larger than for San Marco-2 while the shift of (ji 2 *fi a ) lB 
considerably less. This behavior is what one would exj>ect considering the dif- 


► 


ferences in the orbits. 
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Figures 3 and 4 are graphs of semi-major axis versus time (days from in- 
sertion), for San Marco-2 and Cannonball during an actual lifetime study for the 
two satellites. The circles are those values of a, predicted by the BMW theory 
using only the initial conditions given above, while the crosses indicate those 
arcs of data supplied by orbit improvement routines utilizing tracking or ob- 
servational data to update the orbital elements. The orbital improvements Were 
necessitated by the rapid deterioration of the orbital parameter quality due to 
inaccurate force modeling (particularly air resistance accelerations) in the 
equations of motion. The accepted date of re-entry into the earth's atmosphere 
for San Marco-2 was on October 14, 19G7 at approximately 13 hours. Thus, the 
total lifetime was about 171 days. The BMW program computed a lifetime of 
165 days, for a re-entry on October 8, 1967. Cannonball's re-entry date was 
approximately January 28, 1972, a lifetime of about 173 days. BMW computed 
170 days. In these two cases, the program evaluated the limits of the variational 
equations for values of the eccentric anomaly corresponding to five day intervals. 
This was done so as to allow the BMW program to compare at intermediate 
points, its values for semi-major axis with the observed ones. In actual prac- 
tice, these limits would be evaluated for those values of eccentric anomaly 
corresponding to those regions of an orbit over which a fit of the density variation 
remained fixed. Using a change in perigee height criterion (preselected at 1 km), 
results from Cannonball show that the first region is the first 65 days, the 
second is the next 55 days, and so on, until within the last month of life the regions 
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are only of 5 day durations. As a result, the entire ephemeris is computed 
rapidly. Using the IBM 360/91 electronic digital computer, the complete San 
Marco-2 and Cannonball ephemerides computed at 5 day intervals were exe- 
cuted in less than 18 seconds of computer CPU time. 

These results show that the BMW program was within 4 percent of the 
"true" lifetimes of both satellites, despite the fact that Cannonball had a re- 
latively high eccentricity and low perigee. In addition, it must be pointed out 
that the initial conditions used in the program were obtained from orbit im- 
provement routines other than the Vinti Orbit Determination System, since this 
was the only data available. To be more consistent, and to improve accuracy, 
one should utilize, if possible, only those initial or epoch values obtained from 
the Vinti Orbit Determination System (Reference 9, 10). 

In spite of the startling success with San Marco-2 and Cannonball, improve- 
ments to the present BMW computer program are being considered. Prime 
among these is the incorporation of a dynamic model atmospheric density pro- 
file such as Jacchia's. As is well known a slight miscalculation in the exospheric 
temperature for a static model would drastically alter the resulting density 
profile, and consequently, the computed ephemeris of the satellite. 

The distinct advantage of the BMW method over numerical integration tech- 
niques is that being analytic or closed form, it is not subject to large error 
accumulation due to roundoff and truncation. In addition, an entire ephemeris 
can be obtained in a matter of a few seconds on present generation computers, 
while such a task might be prohibitive with a numerical integration. 



With the atmospheric profile gh en by (5), expression (11) with p = 0 could 
be expressed in terms of Bessel functions; for example, 
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exp (c cos E) cos' 1 EdE = T 0 (c) 


(n = 0) 


= (n = 1) 

= .T 0 (O-ij 1 (c) (n = 2) 

and so forth. The disadvantage here, besides the limited step size over which 
to perform the summation, is in the difficulty of handling intermediate points 
within the limits. A Eessel function approach would appear to be better em- 
ployed in a study of atmospheric density inference from calculations of the 
period decrement. 

King-Hele (Reference 6) has studied the contraction of orbits in a closed 
form manner. Four situations are considered; 

1. Normal e. Phase 1: approx. 0.02<e<0.2 (3 < (ae/H p ) <30 approx.) 

2. Normal e. Phase 2: 0<e<0.02 approx. (0 < (ae/H p )< 3) 

3. Circular orbits; e = 0 (ae/H p = 0). 

4. High eccentricity; e> 0.2 (ae/H p i 30 approx.). 

It is felt that the regions covered by the equations of motion in these four cases 

are also covered by the BMW theory. Even though a static model atmosphere 
profile was used in the calculations, BMW does have the latitude of utilizing a 



density profile such as the Jacchia model. In addition, the reference orbit is 
not an exact ellipse,. but the Vinti orbit. As a result, it is felt that the BMW 
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method has contributed to the program for calculating both accurately and rapidly, 

the orbits of satellites experiencing aerodynamic drag. 

At the present time, the Vinti differential correction algorithm (Reference 10) 

employing a classical weighted least squares technique is being appropriately 
modified to accommodate the augmented force model. The bulk of the imple- 
mentation requires the reformulation of the partial derivatives { 3y ; /3q .(t 0 )} 
to reflect air drag where (y 4 } are tracking observables and {q. (t Q )} are the 
set of "epoch" Vinti orbital elements and an atmospheric parameter iteratively 
estimated by the differential correction technique. It -appears that complete 
analyticity can be retained for these partial derivative expressions, thus merging 
the theoretical developments presented here with the practical application of 
orbit estimation. This work will be forthcoming in a future report. 


ACKNOWLEDGMENTS 

The authors are grateful to Dr. John P. Vinti of the Measurement Systems 
Laboratory 7 , Massachusetts Institute of Technology, Cambridge, Massachusetts, 
and to Dr. A. Deprit of the Dept, of Astronomy, of the University of Cincinnati, 
for valuable discussions which contributed to this report. Also to Mr. E. M. Jones 
and Mr. R. A. Gordon of the Goddard Space Flight Center for their assistance 
with some of the numerical calculations of this report. 



36 


w'. 


REFERENCES 

1. Vinti, John P., "Gaussian Variational Equations for Osculating Elements of 
an Arbitrary Separable Reference Orbit," Celestial Mechanics Journal, I, 
No. 3, pp. 367-375, 1973. 

2. Brouwer, D. and G. Hori, "Theoretical Evaluation of Atmospheric Drag 
Effects in the Motion of an Artificial Satellite," Astronomical J., 6(3, No. 2, 
p. 39, 1961. 

3. Sherrill, J. J., "Development of a Satellite Drag Theory Based on the Vinti 
Formulation," (Ph. D. Dissertation), XJniv. of California, Berkeley, 1966. 

4. Vinti. John P<, "Improvement of the Spheroidal Method for Artificial Satel- 
lites," Astronomical J., 74, No. 1, pp. 25-34, 1969. 

5. Sterne, T. E., "An Introduction to Celestial Mechanics," Interscience 
Publishers, Inc., New York, 1960. 

6. King-Hele, Desmond, "Theory of Satellite Orbits in an Atmosphere," 
Butterworth and Co. (Publishers) LTD., London, 1964. 

7. U. S. Standard Atmosphere Supplements, 1966, (ESSA, NASA, USAF), U. S. 
Government Printing Office, Washington, D. C. 

8. Roberts, C. E., "An Analytic Model for Upper Atmosphere Densities Based 
upon Jacchia's 1970 Models," Celestial Mechanics Journal, 4, pp. 368-377, 
1971. 

9. Bonavito, N. L., "Computational Procedure for V inti's Accurate Reference 
Orbit with Inclusion of the Third Zonal Harmonic," NASA TN D-3562, 
August 1966. 

10. Walden, H. and J. S. Watson, "Differential Corrections Applied to Vinti’s 
Accurate Reference Orbit with Inclusion of the Third Zonal Harmonic," 
NASA TN D-4088, August 1967. ,1. ... .. 



TABLE I 


[ Height 

Static ATM Density 

RHO 

DIFF. 

0.205000000000D 

03 

0.301060802633D 03 

0.301060802633D 03 

0.0 

0.206000000 OOOD 

03 

0.292 561 657391D 03 

0.292581 857391 D 03 

0.113686837722D-12 

0. 207000000 OOOD 

03 

0.284373995477D 03 

Q.284373956992D 03 

0.3S48452348G2D-04 

0,2080000 OOOOOD 

03 

0.276427541284D 03 

0.276427541 2 84D 03 

0.1 13688837722D- 12 

0.2 09C90000000D 

03 

0.268733 -99909D 03 

0.2 607331 23768D 03 

0.7014 11933051D- 04 

0.210000000000D 

03 

0.26 128204 0805D 03 

0.261281 90641 5D 03 

0.134389794596D-03 

0.211000000000D 

03 

0.254 0854821 85D 03 

0.254065482 1 85D 03 

0.213162820728D-13 

0.212000000000D 

03 

0.247075276 16 ID 03 

0.247075072 937D 03 

0.203224049809D-03 

0.213000000000D 

03 

0.2 40303494 560D 03 

0.2 40302 8272 3 8D 03 

0.66732 1 502970D- 03 

0,2 H 0000 OOOOOD 

03 

0.2337425! 53 94D 03 

0.233741317168D 03 

0.11 982261 9084D- 02 

0.21 50000 OOOOOD 

03 

0.227385009955D 03 

0.22738337817 ID 03 

0.1 631 78351 055D- 02 

0.2 160000 OOOOOD 

03. 

0.22 12239304 90D 03 

0.2212^20996430 03 

0.183084742385D-02 

0.2 1 7000000000D 

03 

0.215252488448D 03 

0.215250815844D 03 

0.1682 803 92789D- 02 

0.218000000 OOOD 

03 

0.209464193246D 03 

0.2094 63097141D 03 

0.109610524396D-02 

0.21 90000000 OOD 

03 

0.20385274 155 9D 03 

0. 20385274 1559D 03 

0.248689957516D-13 

0.220000000000D 

03 

,O.1984121O7074D 03 

0.1984173013380 03 

-0.5 1 9426438754D-02 

0.22 1000000000D 

03 

0. 193 13648071 9D 03 

0.193156793127D 03 

-0. 2031240822 86D-01 

0.222000000000D 

03 

0.1880202713191) 03 

0.1 68064945792D 03 

-0.446744726403D-01 

0.22 3 0000 OOOOOD 

03 

0. 18395 8096672D 03 

0.18313S719948D 03 

-0.776232757051D-01 

0.229000000000D 

03 

0. 156242 63352 6D 03 

0.156S64335714D 03 

-0.421702188640D 00 

0.233000000000D 

03 

0.14161 1BS0377D 03 

0.1415C6571597D 03 

0.331878055081D-02 

0.2370000 OOOOOD 

03 

0.12830394S660D 03 

0.1283 03 948660D 03 

0.284217094304D-13 

0.245000000000D 

03 

0. 10564 8314722D 03 

0.105642771638D 03 

0.5543083 15159D-02 

0.250000000000D 

03 

0.937586310738D 02 

0.0374877614 120 02 

0.9654932 80362D-02 

0.270000000000D 

03 

0.590301268805D C2 

0.590233812922D 02 

0.67 4558833608D-02 

0.2 790 000 OOOOOD 

03 

0.48290232351 7D 02 

0.482367020943D 02 

0.530257381655D-03 

0.2 8 5 000000 OOOD 

03 

0.42341 51 1471 9D 02 

0.423458469834D 02 

-0.4135521S2309D-02 

0.2 9 900000 OOOOD 

03 

0.3138216458751) 02 

0.31 4306093484 D 02 

-0.48444760932SD-01 

0.342 OOOOOOOOOD 

03 

0.13472001131 6D 02 

' 0.134712759807D 02 

0.7851 50904 545D-03 

0.410000000000D 

03 

0.398B71539430D 01 

0.408449712729D 01 

-0.95781 7329940D-01 

0.480000000000D 

03 

0. 129181 63S981D 01 

0. 1291 5955431 5D 01 

0.221 398651 137D- 03 

0.550000000090D 

03 

0.449473286944D 00 

0.4534502972641) 00 

-0.397701031921D-02 

0.650000000000D 

03 

0.1J20757S0808D 00 

0.1 12756B96810D 00 

-0.68114800J936D-03 





H.awI 



Figure 2. Variation of a, e, and 
(/Sj + / 5 3 ) for Cannonball 
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Figure 3. Variation of Semi-Major Axis for San Marco-2 
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Figure 4. Variation of Semi-Major Axis for Cannonball 





Figure 5* BMW Computational Algorithm for Drag-Oblateness Interaction 
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APPENDIX 

DEVELOPMENT OF VARIATIONAL EQUATION CONSTANTS 

For k = l f 2, . . 6; the general expression for the change in the Vinti 

orbital elements due to atmospheric drag between the arbitrary limits E*, E* is 
given by: 




P- ^ f* j 4 1 

= m ' + 1 J )+ {t j + i exp fd j + i } dG «<E) 

1=0 


Cl-r-x 

+ / m 2n-j I ’■ {t 2n-j eX P [d 2„-j E]} dG k< E > 

7* 


(A-l) 


♦J 


Ej mod 2*rr 


t* exp Id* E] dG k (E) 


j V t * exp [d* E] dG k (E) 

**E* mod 2 it ' 


Without loss of generality, assume that the fitting process is performed over 
the interval [ 0 ,tt) i.e. 

0 = x 0 < <x 2 < ••• < x n _j < x„ =tt, , 

the coefficients t^ are expressed by 


4 
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\-l~ci [1 +b(L, U)x 2 (cos E l - cos E) 2 ] {l - 1, 2 ••• n) (A-2) 

where is a previously defined fitted parameter in the interval , x^] 

£ - 1, 2, . . . n; x = ae. From symmetry considerations define 


■2nU = 2 ”~*t 

(A-3) 

= (277 - 277 _ X^_ j) 

(A-4) 

L, U) x 2 (cos E^~ cos E) 2 ] l = 1. 2, . 

. . n (A-5) 

d 2n-^ + l = - 

(A-6J 

c 2„.-e + i = H 

(A-7) 


where d^ is the other previously defined, subinterval dependent, fitted parameter. 

The remainder of this appendix shows the procedure utilized to derive the 
values of the non negative integer multiplicity factors m^, m n+ ^ 4 = 1 , 2 , ... n; 
the value f , an integer between 0 , 2 n - 2 ; the value oi , an integer between 
1 , 2 n - 1 ; and (t*. t*) which may be (t^, t^), (t', t^), (t^, tp, or (t', t') depending 
upon whether E* mod 2-.r and E* mod 2 77 are less than or equal to or greater 
t han 77, for arbitrary E*, E* and for different modes of operation of the orbit 
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program. Prior to defining this procedure, consider the following definitions 
concerning (A) modes of operation of the Vinti orbit generator program, (B) 
integration intervals of the variational equations, and (C) Fortran variable 
definitions. 

(A) Mode I 

The Vinti orbit generator program is operating in mode I when the time 
span of the variational equation integration [or implicitly, the integration 
interval (E*,E*)] is sufficiently small so that there is not a complete fitting 
subinterval imbedded between the integration limits. Mode I will be the 
dominant mode for definitive ephemerides or an ephemeris computed to support 
a differential correction. 

Mode n 

The Vinti orbit generator is operating in mode n when the time span of 
variational integration has one or more imbedded fitting subintervals be- 
tween the limits of variational equation integration. This mode will be 
utilized in lifetime studies when large intervals of integration will be per- 
formed. Note that mode n implies the assumption of a valid fit (or x = 
ae constant) over large periods of time than does mode I where new fits will 
be performed as frequently as is required. 

(B) Class I Interval 

A class I interval is an interval such that both boundary points are adjacent 
grid points of the mesh x , . . . , X n . 
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C lass II Interval 

A class n intei’val is an interval such that one boundary point of the interval 
is a point x . of the mesh x 0 , x u . . . , x n while the other boundary point is a 
non mesh point between 


and x jt r 


(C) Fortra n Variable Definitions 
Fortran Variables 
.El 
E2 

NOINTV 


Mathematical Description 


2n 


/ I = 2. 


NOINTV 


COEFF (I , J) J = 1 



c , toe 

1 n 


C n*l t0 C 2n 


J - 2 


T „ NO'NTV 

1 = 2 . — + 


fNOINTV ^ 

.J"* - * ' 

l^NOINTV + 1 


d j to d f 


d n + l t0 d 2n 


J = 3 {1 = 1. NOINTV + 1 


x 0 , Xj. ••• x 


2n 



.ARRAY (I, J) 


r-< 

II 

E* mod 27? 

J = 2 

X 7) 

J = 3 

0 or 1 

J - 4 

**, 

J = 5 

d \ 

J = 1 


J = 2 

E* mod 277 

J = 3 

0 or 1 


J = 4 tj 

w J 5 d* 


1 = 3, NOINTV + 


C 


I , 3, NO INTV , ! 


x o> x r •••. x 2n _, 

X l* X 2’ X 2n 

{ J = 3 mji 

J = 4 c l 

J = 5 d Z 


J = 3 m n+ £ 

I = N0 HH + 2 , NOINTV + 2 j J = 4 c' n+ ^ 

J=5 
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The following describes the subtasks associated with defining the full set of 
integration parameters in (A-l). The detailed logic is presented in the flow- 
chart given in Figure A-l . 

(A) This algorithm will test whether E* and E* are within the same sub- 
interval, will store the proper interval of integration and associated constants 
into ARRAY, and exit. This simplified logic is most useful in Mode I orbit 
generations. 

(g) This algorithm defines the Clas3 EE intervals associated with the arbi- 
trary limits E*, E*, (i.e. [x^ , E* mcd2-7r] and [E* mod 2m, x^] ) and stores them 
and the associated constants into ARRAY. 

© This algorithm checks to see if any of the Class II intervals are near 
77/2 or 3 77/2 in order to avoid underflow difficulties during evaluation of the 
integrals. 

@ To assist the task of deleting unnecessary computations in Mode I gen- 
erations, this algorithm establishes the necessary intervals and associated con- 
stants when Ef and E* are in adjacent subintervals and in the same multiple 
of 2 77 * 

(E) This logic stores the entire set of Class I subintervals and associated 
constants into ARRAY and checks the Class n subintervals for negligible length. 
If such a Class II interval is found, it is el<r.. mated from consideration. An 
arbitrary interval length of 1 x 1C " 7 is the criterion presently used. 
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(T) This algorithm defines m n+ ^ , It is predominantly used in Mode II 
orbit generations. 

(8) This algorithm redefines ARRAY to eliminate Class I intervals with a 
multiplicity factor of zero and Class II intervals whose multiplicity factor lias 
been reset from one to zero by (E) wh9n its length is negligible. 




















